!--------------------------------------------------------
! $Id: x2y2model.shar,v 1.2 2011/05/14 21:34:35 daw Exp $
! Description: do equilibrium MC (one routine for step, one for steps)
! Author:  Murray S. Daw <daw@clemson.edu>
! Created: 14 may 11
!--------------------------------------------------------

module emc

  use prec
  use rannums
!  use global_data
  IMPLICIT NONE

  real(r8) :: dq          ! scale of attempted displacements
  real(r8) :: dp          ! scale of attempted change in speed
  integer  :: seed = -123456789

  real(r8), dimension(2) :: qsc
  real(r8), dimension(2) :: psc
  real(r8) :: rand, mr

  real(r8) :: emccurrent

  integer(i12) :: nattempteddq   ! # of attempted displacements in q
  integer(i12) :: naccepteddq    ! # of accepted displacements in p
  integer(i12) :: nattempteddp   ! # of attempted displacements in q
  integer(i12) :: naccepteddp    ! # of accepted displacements in p
  real(r8) :: ratioaccepteddq  ! ratio of accepted in q 
  real(r8) :: ratioaccepteddp  ! ratio of accepted in p
  real(r8) :: lambda,lambda1  ! anharmonicity
  
  
  private :: dq, dp, seed, qsc, psc, rand, emccurrent
  public :: setupemc, stepemc, manystepsemc, closeemc

contains

  subroutine setupemc( dqin, dpin, qs, ps, seedin, lambdain, &
		      lambda1in,  emccurrentout )
    real(r8), intent(in) :: dqin
    real(r8), intent(in) :: dpin
    real(r8), dimension(2), intent(in) :: qs
    real(r8), dimension(2), intent(in) :: ps
    integer, intent(in)  :: seedin
    real(r8), intent(in) :: lambdain, lambda1in
    real(r8), intent(inout) :: emccurrentout

!    mr = mass

    if ( dqin .lt. 0 ) then
       dq = 0.8 
       print *,' dq set to ',dq
    else
       dq = dqin
    endif
    if ( dpin .lt. 0 ) then
       dp = 2.0 
       print *,' dp set to ',dp
    else
       dp = dpin
    endif

    seed = seedin
    rand = rgauss(seed)

    lambda = lambdain
    lambda1 = lambda1in
    
    emccurrent = 0.5*(qs(1)**2 + qs(2)**2 + ps(1)**2 + ps(2)**2) + &
              0.25*lambda*(qs(1)**2 + qs(2)**2)**2 + &
              1./3.*lambda1*(qs(1)**3-3*qs(1)*qs(2)**2)
    emccurrentout = emccurrent

!    open(unit=50,file='attempts.dat',status='unknown')
!    open(unit=51,file='accepted.dat',status='unknown')

    nattempteddq = 0
    nattempteddp = 0
    naccepteddq = 0
    naccepteddp = 0
    
    return
  end subroutine setupemc

  subroutine stepemc( qs, ps, Qaccept, emccurrentout, jumptype )
    
    real(r8), dimension(2), intent(inout) :: qs
    real(r8), dimension(2), intent(inout) :: ps
!    real(r8), intent(in) :: m
    logical, intent(inout), optional :: Qaccept
    real(r8), intent(inout), optional :: emccurrentout
    integer, intent(inout), optional :: jumptype

! local vars
    integer  :: ipick, jpick
    real(r8) :: ebefore
    real(r8) :: eafter
    real(r8) :: de
    real(r8) :: deltap
    real(r8) :: deltaq
    real(r8) :: demc
    real(r8) :: olddqp, olddqm, newdqp, newdqm
    real(r8) :: oldphi, newphi
    real(r8) :: deerr
    logical :: Qacceptlocal
    integer :: jumptypelocal
    real(r8) :: etry    !  this is for debugging
!    real(r8), dimension(2) :: mass

! pick a member of the chain
    ipick = ceiling(2*ran(seed))
!  ps represent momentum NOT velocity
    qsc = qs
    psc = ps
    rand = ran(seed)  ! decide whether to change q or p
    if ( rand .gt. 0.5 ) then  ! change p
       jumptypelocal = 1
       deltap = dp*rgauss(seed)
       psc(ipick) = psc(ipick) + deltap
	demc = (ps(ipick) + deltap/2. )*deltap
    else        ! change q
       jumptypelocal = 2
       oldphi = 0.5*(qsc(1)**2 + qsc(2)**2) + &
              0.25*lambda*(qsc(1)**2 + qsc(2)**2)**2 + &
              (1./3.)*lambda1*(qsc(1)**3-3*qsc(1)*qsc(2)**2)
       deltaq = dq*rgauss(seed)
       qsc(ipick) = qsc(ipick) + deltaq
       jpick = mod(ipick,2)+1
       newphi = 0.5*(qsc(1)**2 + qsc(2)**2) + &
              0.25*lambda*(qsc(1)**2 + qsc(2)**2)**2 + &
              (1./3.)*lambda1*(qsc(1)**3-3*qsc(1)*qsc(2)**2)
       demc = newphi - oldphi
!       print *, 'change of energy = ', demc
    endif

!    print *, "qs(1) in emc =", qs(1)

    Qacceptlocal = .false.
    if ( (demc .lt. 0.) .or. (ran(seed) .lt. exp(-demc))) Qacceptlocal = .true.

    etry = emccurrent + demc

    if ( jumptypelocal .eq. 1 ) then
       nattempteddp = nattempteddp + 1
       if ( Qacceptlocal ) naccepteddp = naccepteddp + 1
    endif
    if ( jumptypelocal .eq. 2 ) then
       nattempteddq = nattempteddq + 1
       if ( Qacceptlocal ) naccepteddq = naccepteddq + 1
    endif
    
    if ( Qacceptlocal ) then
       qs = qsc
       ps = psc
       emccurrent = emccurrent + demc
    endif

    if ( present( Qaccept ) ) Qaccept = Qacceptlocal
    if ( present( emccurrentout ) ) emccurrentout = emccurrent
    if ( present( jumptype ) ) jumptype = jumptypelocal
    return
  end subroutine stepemc

! do several steps, not collecting info
  subroutine manystepsemc( qs, ps, nsteps )

    real(r8), dimension(:), intent(inout) :: qs
    real(r8), dimension(:), intent(inout) :: ps
    integer, intent(in) :: nsteps

! local vars
    integer :: istep
!

    do istep = 1,nsteps
       call stepemc( qs, ps )
    enddo

    return
  end subroutine manystepsemc
    
  subroutine closeemc

    ratioaccepteddq = float(naccepteddq)/float(nattempteddq)
    ratioaccepteddp = float(naccepteddp)/float(nattempteddp)

    print *,'   ' 
    print *,' MC statistics '
    print 9110
9110 format(1x,'                p-type          q-type    ')
    print 9111,nattempteddp,nattempteddq
9111 format(1x,'attempted:',i10,5x,i10)
    print 9112,naccepteddp,naccepteddq
9112 format(1x,' accepted:',i10,5x,i10)
    print 9113,ratioaccepteddp,ratioaccepteddq
9113 format(1x,'    ratio:',2g15.2)
    print *,'   '
    
    return
  end subroutine closeemc

end module emc


